#include "fd_rng.h"
#include <math.h> /* FIXME: ELIMINATE DEPENDENCE AS THIS IS NOT GUARANTEED BIT LEVEL IDENTICAL BETWEEN PLATFORMS */

float
fd_rng_float_robust( fd_rng_t * rng ) {

  /* Generate a float as though a continuum rand was generated in [1,2]
     and then round to nearest even rounded to the closest exactly
     representable float point value in [1,2].  (This can generate all
     such representable values in that range.) */

  ulong u = (ulong)fd_rng_uint( rng );    /* 32-bit uniform rand */
  ulong m = (u >> 9) + (u & 1UL);         /* [0,2^23] trapezoidal rand */
  float f = 1.f + ((float)m)*FLT_EPSILON; /* Exact */

  /* Now apply an exponent to yield the correct final result */

  return ldexpf( f, -(int)fd_rng_coin_tosses( rng ) );
}

float
fd_rng_float_exp( fd_rng_t * rng ) {
  ulong u = 1UL + (fd_rng_ulong( rng ) >> 1);    /* In (0,2^63] uniform */
  float f = ((float)u) * (1.f/(float)(1UL<<63)); /* In uniform (0,1] with reasonably low quant, exact */
  return -logf( f );                             /* Convert to exp */
}

#if 0
/* USE THIS WITH LOGF AT END REPLACED WITH FAST APPROX REJECTION IF NO
   LIBM.  USE TRANSFORMATION ON CONTINUUM RAND IF WANTING NO QUANT
   ARTIFACTS.  NEARLY EQUIVALENT CAN USE:

        F(x) = 1 - exp(-x) = 1 - exp2(-x*log2(e))
     -> y = x*log2(e), x = y*(1/log2(e))
        F(y) = 1 - exp2(-y)
     -> Y = -log2( U )

   where U is a uniform rand in [0,1].  But U can be generated by:

        U = (2^-E)( 1 + M )

   where E is IID integer rand with probability:

        P(E) = 2^-E for E>0
             = 0    otherwise

   and M is an IID uniform rand in [0,1].  Thus:

     -> Y = E - log2( 1 + M )

   yielding a variant with minimal tail artifacts.  In single precision:

      e = fd_rng_coin_tosses( rng );
      m = fd_rng_uint( rng );        
      y = e - lg1p( m / 2^32 );
      return y*(1/log2(e));

   and for double precision:

      e = fd_rng_coin_tosses( rng )
      m = fd_rng_ulong( rng );
      y = e - lg1p( m / 2^64 );
      return y*(1/log2(e));

   Eliminating floating point rep quant artifacts near the origin is
   much trickier due to the cancellations implied above (but also
   practically much less important ... they are no worse than the
   quantization artifacts of the typical generator). */

float
fd_rng_float_exp( fd_rng_t * rng ) {

  for(;;) {

    /* Generate the 7-bit ziggurat level and a 24-bit trapezoidal rand
       (e.g. "uniform" rand on [0,2^24] with the endpoints 0 and 2^24
       half as likely as the other points). */

    ulong u = (ulong)fd_rng_uint( rng );

    ulong l = (u >> 1) & (zig_level_cnt-1UL); /* uniform level (level_cnt must be a power-of-2 <= 128) */
    ulong m = (u >> 8) + (u & 1UL);           /* ~2^24 bit trapezoidal rand */

    /* Compute a uniform rand as wide as the current ziggurat level
       (will be closed at both endpoints due to trapezoid rounding
       above).  If willing to have an additional table, the
       multiplication by 2^-24 can be eliminated (technically though this
       is exact and can be done just be playing with with the exponent). */

    float x = zig_x[l+1UL]*((1.f/(float)(1<<24))*(float)m);

    /* If this is smaller than the level above this, we can immediately
       accept this rand.  This occurs the vast majority of the time. */

    if( FD_LIKELY( x<=zig_x[l] ) ) break; /* Quick accept */

    /* We can't do a quick acceptance */

    if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {

      /* We are picking a point uniformly in the l-th ziggurat strip. */

      float t = zig_y[l];
      float y = t + (zig_y[l+1UL]-t)*fd_rng_float_c0( rng );
      if( y < expf(-x) ) break; /* Fundamentally ~50/50 unpredictable */

    } else {

      /* We are picking a point uniformly in the ziggurat exponential
         tail.  Note that x - R is an independent uniform rand between 0
         and 1.  In finite precision arithmetic, (x-R)==0 will never
         occur due to the above acceptance check.  Thus, it is safe to
         compute ln (x-R) to generate an exponentially distributed rand
         for the tail here. */

      x = zig_r - logf( x - zig_r );
      break;

    }

    /* Rejected */
  }

  return x;
}
#endif

float
fd_rng_float_norm( fd_rng_t * rng ) {

  /* BEGIN AUTOGENERATED CODE *****************************************/

  static float const zig_x[65] = {
    0.000000000000000000000000000000e+00f /* l= 0 */, 3.455520224239806172095417630130e-01f /* l= 1 */,
    4.621556810233836948395482607799e-01f /* l= 2 */, 5.450460069613645023634852793126e-01f /* l= 3 */,
    6.119536252363746848246563170282e-01f /* l= 4 */, 6.693178808991654190761917686547e-01f /* l= 5 */,
    7.202743599269989995133832427765e-01f /* l= 6 */, 7.666087419009196196129432565591e-01f /* l= 7 */,
    8.094445383821111665470331153482e-01f /* l= 8 */, 8.495396579691822718460002261676e-01f /* l= 9 */,
    8.874325556154150376969927394022e-01f /* l=10 */, 9.235214816642028361095943800319e-01f /* l=11 */,
    9.581106750231991023205452284728e-01f /* l=12 */, 9.914388611066755936521502357017e-01f /* l=13 */,
    1.023697658831656335945772817730e+00f /* l=14 */, 1.055043930218579828852511204307e+00f /* l=15 */,
    1.085608336108546385250818444579e+00f /* l=16 */, 1.115501429269408685757483667977e+00f /* l=17 */,
    1.144818099662757484029536325654e+00f /* l=18 */, 1.173640887904643574190903521082e+00f /* l=19 */,
    1.202042503654549259546786832420e+00f /* l=20 */, 1.230087774542608308986184340039e+00f /* l=21 */,
    1.257835180410570578152240628356e+00f /* l=22 */, 1.285338081362219465062293743962e+00f /* l=23 */,
    1.312645717221077709701594626868e+00f /* l=24 */, 1.339804034974004910067517382100e+00f /* l=25 */,
    1.366856386251642895349189821275e+00f /* l=26 */, 1.393844126728532207214715510357e+00f /* l=27 */,
    1.420807142148611732640484106582e+00f /* l=28 */, 1.447784320603175013198252174540e+00f /* l=29 */,
    1.474813987119431635871430463780e+00f /* l=30 */, 1.501934314168690078094073325765e+00f /* l=31 */,
    1.529183720117960190352200677832e+00f /* l=32 */, 1.556601266765526946100374472426e+00f /* l=33 */,
    1.584227066827416976904988055175e+00f /* l=34 */, 1.612102712541159403396816285348e+00f /* l=35 */,
    1.640271737439234119525742483514e+00f /* l=36 */, 1.668780124881012518812430089898e+00f /* l=37 */,
    1.697676879240338253859027295434e+00f /* l=38 */, 1.727014678920001637717554499041e+00f /* l=39 */,
    1.756850634895277701885175913876e+00f /* l=40 */, 1.787247184704316305928380181900e+00f /* l=41 */,
    1.818273160330322024778176848159e+00f /* l=42 */, 1.850005080182403822074664601072e+00f /* l=43 */,
    1.882528731753153053690014173682e+00f /* l=44 */, 1.915941134586538274016051519588e+00f /* l=45 */,
    1.950353006115250048047894682046e+00f /* l=46 */, 1.985891900706324235165341207665e+00f /* l=47 */,
    2.022706262852767138214066244828e+00f /* l=48 */, 2.060970741899825969985216023161e+00f /* l=49 */,
    2.100893279890740158706921580922e+00f /* l=50 */, 2.142724743933560004960706124599e+00f /* l=51 */,
    2.186772297656567175727290730514e+00f /* l=52 */, 2.233418418575520170057252533624e+00f /* l=53 */,
    2.283148713210882645120378131587e+00f /* l=54 */, 2.336593955757456149034331782666e+00f /* l=55 */,
    2.394596149760118465339014948157e+00f /* l=56 */, 2.458317361057721973189790776182e+00f /* l=57 */,
    2.529429815978553923233249078883e+00f /* l=58 */, 2.610473649022010540245164467166e+00f /* l=59 */,
    2.705599986069530120576243081842e+00f /* l=60 */, 2.822342497323976174845167053107e+00f /* l=61 */,
    2.976823798790561124028714035106e+00f /* l=62 */, 3.215929245508522913137711141118e+00f /* l=63 */,
    3.526881360327788386254538322007e+00f /* l=64 */
  };

  static float const zig_y[65] = {
    1.000000000000000000000000000000e+00f /* l= 0 */, 9.420441848916480126344477619149e-01f /* l= 1 */,
    8.987108451876446194939510037081e-01f /* l= 2 */, 8.619676182560859538047716432718e-01f /* l= 3 */,
    8.292416921465940136417270556191e-01f /* l= 4 */, 7.993205594629470858670725053052e-01f /* l= 5 */,
    7.715162251201959925242870874662e-01f /* l= 6 */, 7.453924046791949711084253327176e-01f /* l= 7 */,
    7.206510565419278137484770940802e-01f /* l= 8 */, 6.970774082324479442845238663651e-01f /* l= 9 */,
    6.745103421549214010562514620695e-01f /* l=10 */, 6.528251409771071645448889397834e-01f /* l=11 */,
    6.319228072029464656281239065549e-01f /* l=12 */, 6.117231258029589507463179287594e-01f /* l=13 */,
    5.921599774953054971587604327077e-01f /* l=14 */, 5.731780673128809060960613119828e-01f /* l=15 */,
    5.547305771308242055075976573164e-01f /* l=16 */, 5.367774409030761823035174384877e-01f /* l=17 */,
    5.192840512302357241190831071975e-01f /* l=18 */, 5.022202719018961842071223367068e-01f /* l=19 */,
    4.855596720803138878681283474581e-01f /* l=20 */, 4.692789240423391189252835808965e-01f /* l=21 */,
    4.533573236341009112014902027177e-01f /* l=22 */, 4.377764041761660086522593704483e-01f /* l=23 */,
    4.225196225029520580892660602812e-01f /* l=24 */, 4.075721013736326630086392874830e-01f /* l=25 */,
    3.929204164392402009802710699526e-01f /* l=26 */, 3.785524188002748156032031823237e-01f /* l=27 */,
    3.644570862756677649099863736115e-01f /* l=28 */, 3.506243980520670366787892163751e-01f /* l=29 */,
    3.370452285453843354473546511940e-01f /* l=30 */, 3.237112571905811243399165438861e-01f /* l=31 */,
    3.106148915554719103920833234156e-01f /* l=32 */, 2.977492017031598617858058342112e-01f /* l=33 */,
    2.851078641441273595351531267017e-01f /* l=34 */, 2.726851140512668792348013879767e-01f /* l=35 */,
    2.604757046803618125588395543213e-01f /* l=36 */, 2.484748731607814430951243836465e-01f /* l=37 */,
    2.366783120090016554051349020882e-01f /* l=38 */, 2.250821458811925040717320800621e-01f /* l=39 */,
    2.136829132292288532124100927656e-01f /* l=40 */, 2.024775526650522966412917846846e-01f /* l=41 */,
    1.914633939793006458976280803608e-01f /* l=42 */, 1.806381539102589138648434843870e-01f /* l=43 */,
    1.699999369289590334732324358735e-01f /* l=44 */, 1.595472415092518358135931233477e-01f /* l=45 */,
    1.492789726065822738132905095343e-01f /* l=46 */, 1.391944614029269361199806984142e-01f /* l=47 */,
    1.292934938280917421830808894390e-01f /* l=48 */, 1.195763500012655360156695570628e-01f /* l=49 */,
    1.100438576497440733127446653439e-01f /* l=50 */, 1.006974639150281902163835620612e-01f /* l=51 */,
    9.153933202201729600920110732631e-02f /* l=52 */, 8.257247254089309598817534793791e-02f /* l=53 */,
    7.380092428122808796957937671479e-02f /* l=54 */, 6.523000887995579990570474762657e-02f /* l=55 */,
    5.686669921543156531844320777935e-02f /* l=56 */, 4.872017206675432561731060171484e-02f /* l=57 */,
    4.080267659191996105732375653419e-02f /* l=58 */, 3.313098485527892895195507111383e-02f /* l=59 */,
    2.572902254561226234667566450942e-02f /* l=60 */, 1.863323273948142502092796546354e-02f /* l=61 */,
    1.190567663419300047354134047123e-02f /* l=62 */, 5.678316641776698130187294149412e-03f /* l=63 */,
    0.000000000000000000000000000000e+00f /* l=64 */
  };

  static ulong const zig_level_cnt = 64;

  static float const zig_r      = 3.215929245508522913137711141118e+00f;
  static float const zig_rcp_r  = 3.109521148192654732252473981369e-01f;
  static float const zig_half_r = 1.607964622754261456568855570559e+00f;
  static float const zig_tail   = 2.852700996027780249580940719056e+00f;

  /* END AUTOGENERATED CODE *******************************************/

  ulong s;
  float x;

  for(;;) {

    /* Generate the sign bit, up to 6-bit ziggurat level and a 24-bit
       trapezoidal rand (e.g. "uniform" rand on [0,2^24] with the
       endpoints 0 and 2^24 half as likely as the other points). */

    ulong u = (ulong)fd_rng_uint( rng ); /* 32-bit rand */

    /**/  s = (u >> 1) & 1UL;                 /* random sign */
    ulong l = (u >> 2) & (zig_level_cnt-1UL); /* uniform level (zig_level_cnt must be power-of-2 <= 64) */
    ulong m = (u >> 8) + (u & 1UL);           /* 24-bit trapezoidal rand */

    /* Compute a uniform rand as wide as the current ziggurat level */

    x = zig_x[l+1UL]*((1.f/16777216.f)*(float)m); /* Guaranteed in [0,x[l+1]] */

    /* If this is <= level above this, we can immediately accept this
       rand.  This occurs the vast majority of the time. */

    if( FD_LIKELY( x<=zig_x[l] ) ) break; /* Quick accept */

    /* We can't do a quick acceptance ... the y rand might matter */

    float y = fd_rng_float_c( rng ); /* Guaranteed in [0,1] */

    if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {

      /* We are picking a point uniformly in the l-th ziggurat strip. */

      y = zig_y[l]*(1.f-y) + zig_y[l+1UL]*y; /* Guaranteed in [y0,y1] */

    } else {

      /* We are picking a point uniformly distributed in the ziggurat
         exponential tail.

         The simplest way to do this is to transform a new uniform
         rand into an exponential distribution that joins smoothly to
         the base of the Ziggurat.

         We theoretically can do slightly better by noting that x-r
         is an independent uniform rand in ( r==x[l], x[l+1] ].  Thus,
         instead of computing a new uniform rand for the tail x and
         transforming it to the desired exponential rand, we could get
         an exponential rand via:

           x' = r - (1/r) ln ( (x-r) / (x[l+1]-r) )

         which simplifies to:

           x = (r + (1/r) ln( x[l+1]-r) ) - (1/r) ln ( x-r )

         The first term is a compile time constant.  The second term
         doesn't require a new rand but is otherwise the same
         transformation cost.  x-r is guaranteed greater than 0 given
         the quick accept test above so the logf will not blow up. */

      x  = zig_tail - zig_rcp_r*logf( x - zig_r );
      y *= expf( zig_r*(zig_half_r - x) );

    }

    /* FIXME: COULD GET RID OF 0.5 MULTIPLICATION BY USING PDF
       expf(-x^2), ADJUST TAIL COMP AND SCALING TO UNIT NORM AS PART OF
       THE SIGN INSERTION. */

    if( y < expf((-0.5f)*(x*x)) ) break; /* Slow accept, approx 50/50 unpredictable */

    /* Rejected, try again */
  }

  return fd_float_if( (int)s, -x, x ); /* FIXME: LOAD AND MUL BY +/-1?  LOAD AND COPYSIGNF BY +/-1? */
}

#if FD_HAS_DOUBLE /* See float variants for details how these work */

double
fd_rng_double_robust( fd_rng_t * rng ) {
  ulong  u = fd_rng_ulong( rng );
  ulong  m = (u >> 12) + (u & 1UL);
  double d = 1. + ((double)m)*DBL_EPSILON;
  return ldexp( d, -(int)fd_rng_coin_tosses( rng ) );
}

double
fd_rng_double_exp( fd_rng_t * rng ) {
  ulong  u = 1UL + (fd_rng_ulong( rng ) >> 1);
  double d = ((double)u) * (1./(double)(1UL<<63));
  return -log( d );
}

double
fd_rng_double_norm( fd_rng_t * rng ) {

  /* BEGIN AUTOGENERATED CODE *****************************************/

  static double const zig_x[65] = {
    0.000000000000000000000000000000e+00, 3.455520224239806172095417630130e-01,
    4.621556810233836948395482607799e-01, 5.450460069613645023634852793126e-01,
    6.119536252363746848246563170282e-01, 6.693178808991654190761917686547e-01,
    7.202743599269989995133832427765e-01, 7.666087419009196196129432565591e-01,
    8.094445383821111665470331153482e-01, 8.495396579691822718460002261676e-01,
    8.874325556154150376969927394022e-01, 9.235214816642028361095943800319e-01,
    9.581106750231991023205452284728e-01, 9.914388611066755936521502357017e-01,
    1.023697658831656335945772817730e+00, 1.055043930218579828852511204307e+00,
    1.085608336108546385250818444579e+00, 1.115501429269408685757483667977e+00,
    1.144818099662757484029536325654e+00, 1.173640887904643574190903521082e+00,
    1.202042503654549259546786832420e+00, 1.230087774542608308986184340039e+00,
    1.257835180410570578152240628356e+00, 1.285338081362219465062293743962e+00,
    1.312645717221077709701594626868e+00, 1.339804034974004910067517382100e+00,
    1.366856386251642895349189821275e+00, 1.393844126728532207214715510357e+00,
    1.420807142148611732640484106582e+00, 1.447784320603175013198252174540e+00,
    1.474813987119431635871430463780e+00, 1.501934314168690078094073325765e+00,
    1.529183720117960190352200677832e+00, 1.556601266765526946100374472426e+00,
    1.584227066827416976904988055175e+00, 1.612102712541159403396816285348e+00,
    1.640271737439234119525742483514e+00, 1.668780124881012518812430089898e+00,
    1.697676879240338253859027295434e+00, 1.727014678920001637717554499041e+00,
    1.756850634895277701885175913876e+00, 1.787247184704316305928380181900e+00,
    1.818273160330322024778176848159e+00, 1.850005080182403822074664601072e+00,
    1.882528731753153053690014173682e+00, 1.915941134586538274016051519588e+00,
    1.950353006115250048047894682046e+00, 1.985891900706324235165341207665e+00,
    2.022706262852767138214066244828e+00, 2.060970741899825969985216023161e+00,
    2.100893279890740158706921580922e+00, 2.142724743933560004960706124599e+00,
    2.186772297656567175727290730514e+00, 2.233418418575520170057252533624e+00,
    2.283148713210882645120378131587e+00, 2.336593955757456149034331782666e+00,
    2.394596149760118465339014948157e+00, 2.458317361057721973189790776182e+00,
    2.529429815978553923233249078883e+00, 2.610473649022010540245164467166e+00,
    2.705599986069530120576243081842e+00, 2.822342497323976174845167053107e+00,
    2.976823798790561124028714035106e+00, 3.215929245508522913137711141118e+00,
    3.526881360327788386254538322007e+00
  };

  static double const zig_y[65] = {
    1.000000000000000000000000000000e+00, 9.420441848916480126344477619149e-01,
    8.987108451876446194939510037081e-01, 8.619676182560859538047716432718e-01,
    8.292416921465940136417270556191e-01, 7.993205594629470858670725053052e-01,
    7.715162251201959925242870874662e-01, 7.453924046791949711084253327176e-01,
    7.206510565419278137484770940802e-01, 6.970774082324479442845238663651e-01,
    6.745103421549214010562514620695e-01, 6.528251409771071645448889397834e-01,
    6.319228072029464656281239065549e-01, 6.117231258029589507463179287594e-01,
    5.921599774953054971587604327077e-01, 5.731780673128809060960613119828e-01,
    5.547305771308242055075976573164e-01, 5.367774409030761823035174384877e-01,
    5.192840512302357241190831071975e-01, 5.022202719018961842071223367068e-01,
    4.855596720803138878681283474581e-01, 4.692789240423391189252835808965e-01,
    4.533573236341009112014902027177e-01, 4.377764041761660086522593704483e-01,
    4.225196225029520580892660602812e-01, 4.075721013736326630086392874830e-01,
    3.929204164392402009802710699526e-01, 3.785524188002748156032031823237e-01,
    3.644570862756677649099863736115e-01, 3.506243980520670366787892163751e-01,
    3.370452285453843354473546511940e-01, 3.237112571905811243399165438861e-01,
    3.106148915554719103920833234156e-01, 2.977492017031598617858058342112e-01,
    2.851078641441273595351531267017e-01, 2.726851140512668792348013879767e-01,
    2.604757046803618125588395543213e-01, 2.484748731607814430951243836465e-01,
    2.366783120090016554051349020882e-01, 2.250821458811925040717320800621e-01,
    2.136829132292288532124100927656e-01, 2.024775526650522966412917846846e-01,
    1.914633939793006458976280803608e-01, 1.806381539102589138648434843870e-01,
    1.699999369289590334732324358735e-01, 1.595472415092518358135931233477e-01,
    1.492789726065822738132905095343e-01, 1.391944614029269361199806984142e-01,
    1.292934938280917421830808894390e-01, 1.195763500012655360156695570628e-01,
    1.100438576497440733127446653439e-01, 1.006974639150281902163835620612e-01,
    9.153933202201729600920110732631e-02, 8.257247254089309598817534793791e-02,
    7.380092428122808796957937671479e-02, 6.523000887995579990570474762657e-02,
    5.686669921543156531844320777935e-02, 4.872017206675432561731060171484e-02,
    4.080267659191996105732375653419e-02, 3.313098485527892895195507111383e-02,
    2.572902254561226234667566450942e-02, 1.863323273948142502092796546354e-02,
    1.190567663419300047354134047123e-02, 5.678316641776698130187294149412e-03,
    0.000000000000000000000000000000e+00
  };

  static ulong const zig_level_cnt = 64;

  static double const zig_r      = 3.215929245508522913137711141118e+00;
  static double const zig_rcp_r  = 3.109521148192654732252473981369e-01;
  static double const zig_half_r = 1.607964622754261456568855570559e+00;
  static double const zig_tail   = 2.852700996027780249580940719056e+00;

  /* END AUTOGENERATED CODE *******************************************/

  ulong  s;
  double x;

  for(;;) {

    ulong u = fd_rng_ulong( rng );             /* 64-bit rand */

    /**/  s = (u >>  1) & 1UL;
    ulong l = (u >>  2) & (zig_level_cnt-1UL); /* uniform level (zig_level_cnt must be power-of-2 <= 512) */
    ulong m = (u >> 11) + (u & 1UL);           /* 53-bit trapezoidal rand */

    x = zig_x[l+1UL]*((1./9007199254740992.)*(double)m);

    if( FD_LIKELY( x<=zig_x[l] ) ) break;

    double y = fd_rng_double_c( rng );

    if( FD_LIKELY( l<(zig_level_cnt-1UL) ) ) {
      y = zig_y[l]*(1.-y) + zig_y[l+1UL]*y;
    } else {
      x  = zig_tail - zig_rcp_r*log( x - zig_r );
      y *= exp( zig_r*(zig_half_r - x) );
    }

    if( y < exp((-0.5)*(x*x)) ) break;
  }

  return fd_double_if( (int)s, -x, x );
}

#endif

